Prognostic value and immune landscapes of immunogenic cell death-associated lncRNAs in lung adenocarcinoma

Immunogenic cell death (ICD) has been demonstrated to activate T cells to kill tumor cells, which is closely related to tumor development, and long noncoding RNAs (lncRNAs) are also involved. However, it is not known whether ICD-related lncRNAs are associated with the development of lung adenocarcinoma (LUAD). We downloaded ICD-related genes from GeneCards and the transcriptome statistics of LUAD patients from The Cancer Genome Atlas (TCGA) and subsequently developed and verified a predictive model. A successful model was used together with other clinical features to construct a nomogram for predicting patient survival. To further study the mechanism of tumor action and to guide therapy, we performed enrichment analysis, tumor microenvironment analysis, somatic mutation analysis, drug sensitivity analysis and real-time quantitative polymerase chain reaction (RT-qPCR) analysis. Nine ICD-related lncRNAs with significant prognostic relevance were selected for model construction. Survival analysis demonstrated that overall survival was substantially shorter in the high-risk group than in the low-risk group (P < 0.001). This model was predictive of prognosis across all clinical subgroups. Cox regression analysis further supported the independent prediction ability of the model. Ultimately, a nomogram depending on stage and risk score was created and showed a better predictive performance than the nomogram without the risk score. Through enrichment analysis, the enriched pathways in the high-risk group were found to be primarily associated with metabolism and DNA replication. Tumor microenvironment analysis suggested that the immune cell concentration was lower in the high-risk group. Somatic mutation analysis revealed that the high-risk group contained more tumor mutations (P = 0.00018). Tumor immune dysfunction and exclusion scores exhibited greater sensitivity to immunotherapy in the high-risk group (P < 0.001). Drug sensitivity analysis suggested that the predictive model can also be applied to the choice of chemotherapy drugs. RT-qPCR analysis also validated the accuracy of the constructed model based on nine ICD-related lncRNAs. The prognostic model constructed based on the nine ICD-related lncRNAs showed good application value in assessing prognosis and guiding clinical therapy.


IRG
Immunogenic-cell-death related gene KEGG Kyoto Encyclopedia of Genes and Genomes K-M survival analysis Kaplan Accounting for 21% of all cancer deaths globally, lung cancer (LC) has become the primary cause of mortality 1 .
The most common subtype of LC is lung adenocarcinoma (LUAD), and the morbidity and mortality of LUAD continue to increase yearly 2 .Tumor node metastasis (TNM) staging is frequently used to forecast clinical outcomes, but the predictive effect is still unsatisfactory 3,4 .Thus, it is imperative to build a better assessment measure for predicting patient survival and guiding LUAD treatment.In recent years, the approach of constructing a predictive model through a combination of several biomarkers and using it to assess tumor patient prognosis has been widely used [5][6][7] .By activating T cells to produce direct impacts on tumor cell killing and antitumor immune responses, immunogenic cell death (ICD) is an example of regulated cell death that can regulate the growth of tumors 8 .Long noncoding RNAs (lncRNAs) can regulate tumor development by affecting tumor cell metabolism and certain oncogenic or carcinogenic factors; for example, cancer-associated fibroblast-specific lncRNA (LINC01614) enhances glutamine uptake in LUAD, thereby promoting cancer cell growth 9 .A variety of models constructed with lncRNAs to predict the prognosis and treatment options for various cancer types are now available and have also shown good prognostic value 5,6,10 .ICD-related lncRNA models have been developed to forecast the development of high-grade gliomas (HGGs), head and neck squamous cell carcinoma (HNSC) and gastric cancer (GC) [11][12][13] , but there remains a dearth of ICD-related biomarkers for evaluating LUAD prognosis.
In this study, we first used Pearson's analysis to derive ICD-associated lncRNAs that play a role in LUAD, followed by differential analysis, univariate and multivariate cox analysis, and lasso regression to finally screen nine lncRNAs to construct a prognostic model.After the successful verification of this model, we created a nomogram to estimate the survival time of patients with LUAD.Subsequently, we explored the likely mechanisms by which ICD-associated lncRNAs act in LUAD and provided ideas for options for clinical treatment.

Statistics source
We downloaded the transcriptome statistics for LUAD from The Cancer Genome Atlas (TCGA) (https:// portal.gdc.cancer.gov/, until October 29th, 2022), consisting of 59 normal samples and 526 tumor samples (585 samples), including FPKM data and count data.The count data were log2-transformed using the "limma" package 14 .Additionally, for further study, the survival information, clinical information, and simple nucleotide variation (SNP) information for LUAD were retrieved from TCGA.The term "immunogenic cell death" was searched in GeneCards (https:// www.genec ards.org/, until October 29th, 2022), and we selected 138 ICD-related genes (IRGs) based on correlation scores > 35.All statistical analyses were carried out in accordance with R4.2.1.

Investigation of ICD-related lncRNAs with differential expression
Our method is displayed in a flowchart (Fig. 1).To identify ICD-related lncRNAs that are differentially expressed, we isolated lncRNAs expressed in LUAD for Pearson analysis with IRGs, selected ICD-related lncRNAs based on P < 0.05 and |correlation coefficient| > 0.3, and then performed difference analysis with the criteria |log2Fold-Change| > 1 and P < 0.05.

Constructing the model
All cases with normal or no survival details were removed, leaving a total of 477 tumor samples.At a random ratio of 1:1, the 477 patients were divided into two groups, with 239 samples in the training subgroup and 238 samples in the test subgroup.The training cohort was used for model construction.To identify ICD-associated lncRNAs that play a major role in LUAD, we carried out univariate Cox analysis to identify ICD-related lncRNAs linked to prognosis (P < 0.05, HR > 1), followed by LASSO analysis to avoid overfitting.Then, we performed multivariate

Verifying the model's feasibility
The test cohort and entire cohort were used for model validation, and employing the RiskScore algorithm, the risk score for each individual was computed.Then, the median of each cohort was used to classify each cohort into high-and low-risk groups.All three cohorts underwent Kaplan-Meier (K-M) analysis 15 , receiver operating characteristic (ROC) analysis 16 , risk mapping, heatmapping to visualize expression, and principal component analysis (PCA) 17 to assess the predictive effect and determine whether there was a discernible distinction in life expectancy between the two risk groups.Stratified analysis can help us understand whether this prognostic model is applicable to the entire population.To confirm whether the constructed model was an independent predictor, we used age, gender, stage, T stage, N stage, and M stage together with the risk score and conducted univariate and multivariate Cox analyses 18 .

Structuring the nomogram
We investigated the link between the risk score and other clinical factors.To identify whether the risk score has the potential to outperform other clinical indicators as an independent predictor, the predictive effect was

Evaluation of the tumor immune microenvironment
The stromal cell and immune cell concentrations in the tumors were compared between the two risk groups using the ESTIMATE method to determine whether there were any variations.Overall, 16 immune cell types and 13 immune-associated regulatory mechanisms were included in the 29 immunological markers measured and compared by single sample gene set enrichment analysis (ssGSEA) in the two risk groups using the "GSVA" R package 22 .Finally, we constructed an immune-related heatmap.In timer 2.0 (http:// timer.cistr ome.org/, until October 29th, 2022), the proportion of immune cells was determined using seven different methods: TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, XCELL, EPIC, and MCPCOUNTER.Additionally, scatter plots were generated using the CIBERSORT method (P < 0.05) to display the interactions between immune cells and the risk score 23 .

Somatic cell mutation analysis
Simple nucleotide variation (SNP) data were downloaded from TCGA for LUAD.We utilized the "maftools" package to determine the mutation status of individuals and to estimate the tumor mutation burden (TMB) score for every individual.We then compared the TMB between the two risk groups 24 .

Clinical pharmacotherapy
Tumor immune dysfunction and exclusion (TIDE) (http:// tide.dfci.harva rd.edu/, until October 29th, 2022) is a method for simulating immunological escape of tumors and predicts sensitivity to PD1 drugs and CTLA4 drugs, and the probability of immune evasion is positively associated with the TIDE score 25 .We applied the "oncoPredict" R package to predict the susceptibility of LUAD patients to 198 chemotherapeutic agents and to search for the most effective chemotherapeutic drugs for the two risk groups.We defined drugs with P < 0.05 as sensitive drugs 26 .

Real-time quantitative polymerase chain reaction (RT-qPCR) for human lung adenocarcinoma
At the Second Affiliated Hospital of Nanchang University, we collected 18 lung adenocarcinoma tissues from patients after surgical resection.The samples were quickly frozen and kept in liquid nitrogen at − 80 °C.Informed consent was obtained from all participants, and the study was authorized by the ethics committee of the Second Affiliated Hospital of Nanchang University.
All samples were divided into a high-risk group and a low-risk group.According to the manufacturer's recommendations, RNA was obtained from LUAD tissues using TRlzol reagent (Life Technologies CA, USA) and then randomly assigned for RT-qPCR analysis.This was followed by reverse transcription using the SureScript First-Strand cDNA Synthesis kit (GeneCopropol, Guangzhou, China) at 45 °C for 1 h.The Applied Biosystems 7500 Fast Real-Time PCR System and BlazeTaq SYBR Green qPCR Master Mix (GeneCopropol, Guangzhou, China) were used to complete the RT-qPCR analysis (Applied Biosystems).To determine the level of RNA expression in each sample, we employed 2 ΔΔCt values.
In the Human Protein Atlas (HPA) (https:// www.prote inatl as.org/, until January 20th, 2023), we analyzed the variations in protein expression of selected ICD-related genes in normal tissues and LUAD tissues, which provides some evidence for the role of ICDs in LUAD.

Ethics approval and consent
The current study investigated the publicly available data, and was also approved by the ethics committee of the Second Affiliated Hospital of Nanchang University (Nanchang, China).All methods were carried out in accordance with the Declaration of Helsinki.

Informed consent
Written informed consent was obtained from all participants included in the study.

Development of a prognostic model and internal validation
We constructed a prognostic model using 9 filtered ICD-related lncRNAs, and risk scores were determined for every sample in the three cohorts.According to K-M analysis, in all three cohorts, the low-risk group had a longer life expectancy than the high-risk group (P < 0.001) (Fig. 3A-C).The reliability of estimating survival rates for 1 year, 3 years, and 5 years was evaluated using the AUC values in time-ROC curves (training: 0.741, 0.783, 0.905; testing: 0.656, 0.634, 0.564; entire: 0.691, 0.702, 0.706) (Fig. 3D-F).Heatmap illustrating the expression of nine ICD-related lncRNAs modeled among the three cohorts (Fig. 3G-I).People with the disease in the low-risk group had a better chance of surviving than those in the high-risk group, according to risk plots for the three cohorts (Fig. S1A-F).PCA was performed with all genes, all lncRNAs, all ICD-related lncRNAs, and risk ICD-related lncRNAs for model construction, in which only risk ICD-related lncRNAs could separate the two risk groups better (Fig. S1G-J).A clinical correlation heatmap revealed that the expression landscape of the nine risk ICD-related lncRNAs modeled was not related to other clinical features (Fig. S2).In addition, K-M analysis showed that the constructed model was fit to predict the prognosis for all LUAD patients (age ≤ 65, age > 65; female, male; stage I-II, stage III-IV; T1-T2, T3-T4; N0, N1-3; M0, M1) (Fig. S3).S4).
In contrast to age and gender, we found that M stage (P = 0.0021), T stage (P = 0.01), and N stage (P = 0.0011) were considerably linked with the risk score (Fig. 4C-H).Then, we drew ROC curves and obtained AUC values (risk score = 0.727, stage = 0.702, gender = 0.544, age = 0.486), revealing that the risk score showed superior predictive value over other clinical indicators (Fig. 5A).A nomogram was constructed on the basis of stage and risk, and different patients could be scored by this system to allow prediction of rates of survival at 1 year, 3 years, and 5 years (Fig. 5B).Then, we chose a patient for successful validation (Fig. S4).Calibration curves with risk showed that the prediction curve highly overlapped with the best prediction line of 45° (C-index = 0.7291493), which had a better fit than the nomogram without risk (C-index = 0.6799758) (Fig. 5C,D).DCA indicated that the model with the risk score had a better clinical benefit than the model without the risk score (Fig. 5E).

Functional analysis
We performed enrichment analysis using five different pathway libraries in GSEA, and the results were similar.The high-risk group showed particular upregulation of pathways that are engaged in DNA replication and metabolism.For example, in the KEGG database, cellular pathways such as DNA replication, cell cycle and citrate cycle TCA cycle were significantly enriched in the high-risk group.In the low-risk group, the pathways that were engaged in immunity and inflammation were particularly upregulated.For example, in the WP database, cellular pathways such as the T-cell receptor signaling pathway, T-cell activation sarscov2 and Th17 cell differentiation pathway were significantly enriched in the low-risk group (Fig. 6A-F).

Tumor immune microenvironment analysis
High-risk group members presented a relatively lower average stromal score (P = 0.0081), a lower immune score (P < 0.001), and a lower ESTIMATE score (P < 0.001) than low-risk group members (Fig. 7A-C).The low-risk group presented a state of moderate activity and had more immune cells and immune functional pathways (Fig. 7E,F).Tumor purity is the proportion of tumor cells to all cells in the sample.A lower estimate score indicates a lower proportion of stromal and immune cells in the tumor and, conversely, a higher proportion of tumor cells, and in the study, we can conclude that the high-risk group suffered from high tumor purity and low immune-related marker expression (Fig. 7D).Higher concentrations of most immune cells corresponded to lower risk scores based on seven different immune algorithms, in which the CIBERSORT algorithm resulted in P < 0.05 immune cells (Fig. S5A).The majority of eosinophils, activated memory CD4 T cells, M0 macrophages and M2 macrophages had a higher proportion in the high-risk group.The majority of memory B cells, plasma cells, resting memory CD4 T cells and regulatory Tregs had a higher proportion in the low-risk group (Fig. S5B-I).

Somatic mutation analysis
We mapped the first 20 genes with the most frequent mutations in the two risk groups, and the most commonly altered gene in both groups was TP53 (Fig. 8A,B).We found that the top 5 genes with the highest mutation rates in the two risk groups were TP53, TTN, MUC16, CSMD3, and RYR2.Mutations in these five genes may be closely related to LUAD progression.In comparison to the low-risk group, the TMB was substantially higher in the high-risk group (P = 0.00018), indicating that these patients may experience a worse result and respond more favorably to immunotherapy (Fig. 8C).

Immunotherapy and chemotherapy
The TIDE results showed notably lower scores in the high-risk group than in the low-risk group (P < 0.001), suggesting that the high-risk group was more amenable to immunological treatment (Fig. 8D).Eleven anticancer medicines were effective in the high-risk group, primarily by blocking the EGFR and ERK-MAPK signaling pathways (Table S5).The low-risk group was susceptible to 116 anticancer agents, mainly by inhibiting the PI3K/ mTOR signaling pathway and DNA replication (Table S6).In both groups, there were no differences in sensitivity to 68 chemotherapeutic agents (Table S7).

Biological validation
The HPA database was used to compare the protein expression of ICD-related genes in normal tissues and LUAD tissues by cellular immunohistochemical staining.The protein expression of the genes is shown in Fig. S6A.The primer sequences of nine ICD-related lncRNAs modeled for RT-qPCR are shown in Table S8.The RT-qPCR results showed that the nine ICD-related lncRNAs used to construct the model were differentially expressed in the two risk groups, which was consistent with our study (Fig. S6B).

Discussion
One of the principal reasons for human mortality is cancer 27 .More than 350 deaths are caused by LC daily.LUAD is the main type of LC, and its incidence remains high 1 .Traditional TNM staging still has significant limitations for assessing patient prognosis and selecting treatment options 4 .New predictive tools need to be discovered 28 .
Selective induction of cancer cell death is the most effective way to fight cancer 29 .When tumor cells are stimulated externally, the conversion of nonimmunogenicity to immunogenicity mediates the induction of ICD by the body to generate an antitumor immune response 8 .The high-risk group showed active pathways that are associated with DNA replication and metabolism.The immune system of patients in the high-risk group remained suppressed.The comparatively high TMB and low TIDE scores in the high-risk group suggest that these patients might be more susceptible to anti-PD1 and anti-CTLA4 treatments.The two risk groups had dissimilar susceptibilities to distinct chemotherapeutic agents.The construction of prognostic models using ICD-related lncRNAs and the prediction of cancer patient prognosis have been employed for a variety of tumor types, including HGGs, HNSC and GC [11][12][13] .Our findings were consistent with these previous studies, with the high-risk group experiencing a poorer outcome and the nomogram used to assess patient prognosis having a better predictive value, suggesting that prognostic features constructed based on ICD-related lncRNAs are a good predictive tool for assessing patient prognosis.Therefore, to construct a similar prognostic model, we identified nine ICD-related lncRNAs via Pearson analysis, univariate Cox analysis, LASSO regression, and multivariate Cox analysis.To accurately forecast the survival rate, we created a nomogram that was successfully validated.Our model has great scientific validity, and two ICD-related lncRNAs affecting tumor mechanisms were investigated.Signal transducer and activator of transcription is negatively regulated when LINC00908 is knocked down, encouraging the growth of tumors 30 .Based on our constructed model, the high-risk group had low LINC00908 expression and a worse prognosis, consistent with previous findings.KIAA0125 overexpression has also been shown to inhibit tumor cell proliferation, metastasis and infiltration through Wnt/β-linked protein signaling 31 , which is consistent with our study, where high KIAA0125 expression in the low-risk group inhibited tumor proliferation and metastasis, leading to a better prognosis.Subsequently, we explored the molecular mechanisms involved in the progression of tumors by enrichment analysis, and to understand the efficacy of immunotherapy in LUAD patients, we performed tumor microenvironmental analysis, tumor mutation burden, and TIDE.Finally, we predicted the sensitivity of LUAD patients to various chemotherapeutic agents, which provides a basis for the clinical search for appropriate drugs.
GSEA was conducted in the two risk groups to identify potential pathways, and immunological and inflammatory-associated pathways were notably enriched in the low-risk group, which indicated a close correlation between the development of LUAD and immunity.A large percentage of neoplastic cells may be killed when ICD is activated in the low-risk group, improving prognosis 8 .Pathways related to metabolism and DNA replication were more abundant in the high-risk group, among which the pentose phosphate pathway (PPP) not only synthesizes pentose phosphate to provide nucleic acids as a raw material for cancer cell proliferation but also generates nicotinamide adenine dinucleotide phosphate (NADPH), which is necessary for redox reactions  in the metabolism of various substances.Activation of the PPP allows cancer cells to metabolize nutrients and actively proliferate, promoting tumor progression.Some tumor suppressors, such as TP53, are also associated with PPP regulation, and in this study, we astonishingly observed that TP53 showed the highest mutation level in the high-risk group.Moreover, TP53 deletion has been shown to reduce the ability of TP53 to prevent cancer cells from absorbing glucose, contributing to overactivation of the PPP and resulting in a poor outcome 32 .The high-risk group is immunosuppressed, and the immunosuppressive state of most tumors can significantly limit ICD-driven immunity to clear tumor cells.In addition, poor immune cell infiltration, including antigenpresenting cells and their precursor cells, means that dead tumor cells are less likely to be effectively processed and drive ICD, which leads to a poor prognosis 8 .A higher TMB indicates a higher possibility of neoantigens and a higher immune reaction rate 33 .In the high-risk group, the TMB was higher and the TIDE scores were lower.According to these findings, the high-risk group seems to be more responsive to immunological therapies, including anti-PD1 and anti-CTLA4 agents.In contrast, the TIDE scores remained higher in the low-risk group, indicating a higher likelihood of immunological escape.Patients with early-stage lung cancer rely mainly on surgical resection, and in patients with advanced lung adenocarcinoma, EGFR inhibitors usually have good benefits 34 .In our study, the high-risk group was found to be sensitive to EGFR inhibitors, such as erlotinib, lapatinib, and gefitinib, which are already being used as the principal treatments for LC.
We present the following innovative ideas in our study.First, we are the first to develop a prognostic model using ICD-related lncRNAs to predict the outcome of patients with LUAD.Second, our risk model has a better prediction effect than alternative clinical indicators.Moreover, we compared a nomogram constructed using only the traditional stage to a nomogram that included the risk model and found that the latter had a better prediction capacity.Third, we found that Lu et al. developed a pattern of necroptosis-related lncRNAs to predict overall survival in LUAD 35 .We performed somatic mutation analysis, immunotherapy analysis, and drug sensitivity analysis, which they did not perform.Furthermore, the AUC value of their constructed model (0.723) was smaller than that of ours (0.727).Therefore, our model is more useful in assessing the prognosis of LUAD patients and selecting treatments.Although our results have good scientific validity, there are still shortcomings.First, we only chose the TCGA database for the analysis, without external validation.Second, the specific mechanisms of action of the ICD-related lncRNAs in LUAD investigated in our study are still not clear and need to be explored in subsequent experiments.The advancement of interaction prediction research in various fields of computational biology would provide valuable insights into genetic markers and related diseases, such as the NDALMA model and GCNCRF model for predicting interactions between lncRNAs and miRNAs 36,37 , the GFPA model and the scAAGA model for processing single-cell data 38,39 , the MDA-AENMF model and the GCNAT model for forecasting potential relationships between metabolites and disease 40,41 , and the DMFGAM model for developing related drugs 42 .In follow-up studies, we can make predictions by developing or using similar models to investigate the mechanism of action of lncRNAs in LUAD and the clinical development of related drugs.Third, this thesis only considered the role of ICD in LUAD, but there are many other types of cell death within the cell.Similar to what Li et al. found, whether our screened lncRNAs also bind to different substances (e.g., miRNAs, proteins) to drive droplet assembly and mediate different types of cell death modes needs to be explored 43,44 .This may help us to have a more in-depth understanding of disease mechanisms and therapeutic targets.

Conclusion
Based on 9 ICD-related lncRNAs, we developed a predictive model and a nomogram that has great value in assessing prognosis and directing clinical therapy in LUAD patients.The presence of DNA replication and metabolism-related pathways and an immunosuppressed state in the high-risk group likely results in a poorer outcome.The risk model is also valuable in guiding the choice of antineoplastic agents for LUAD patients.The RT-qPCR results also further confirmed the accuracy of this model.Due to the above shortcomings, this study still requires basic research to further confirm the specific mechanisms, and the selection of clinically relevant drugs must be verified through clinical practice.

Figure 2 .
Figure 2. Screening of ICD-related lncRNAs for model construction.(A) The volcano plot of differential analysis; (B) Correspondence between IRG and the nine ICD-related lncRNAs; (C,D) 14 ICD-related lncRNAs were identified by lasso regression; (E) Correlation of nine ICD-related lncRNAs with prognosis in the construct model.

Figure 4 .Figure 5 .
Figure 4. Verify that risk score is an independent predictive indicator.(A) Univariate cox regression analysis; (B) Multivariate cox regression analysis; (C-H) Correlation of different clinical characteristics with risk score, including age, gender, stage, T, N, M.

Figure 7 .
Figure 7. Tumor immune microenvironment analysis.(A-C) Comparing the differences in stromal score, immune score, and estimate score between two risk groups; (D) Differential infiltration of tumor immune microenvironment in two risk groups; (E,F) Boxplot displaying the expression of 16 immune cells and 13 immune function pathways in two risk groups based on ssGSEA algorithm.

Figure 8 .
Figure 8. Somatic mutation analysis and immunotherapy.(A,B) Top 20 genes with the highest mutation frequency in tumor cells of the two risk groups; (C) Comparison of TMB between two risk groups; (D) Comparison of TIDE scores between two risk groups.

Table 2 .
The classification of the high-risk group and low-risk group was performed utilizing the middle value for each cohort, as shown for the entire cohort grouping in

Table 1 .
Clinical information of 477 LUAD samples in the TCGA database.LUAD lung adenocarcinoma, TCGA The Cancer Genome Atlas, T tumor, N node, M metastasis.

Table 2 .
Clinical features of LUAD patients in two risk groups.LUAD lung adenocarcinoma.